/*	Interactive Gravitational Simulator
	Project: IGS Core
	File: Gravity.cpp
	Author: Mike Bantegui <mbante2@gmail.com>, Hofstra University 
	Copyright (C) 2012 - 2013 Mike Bantegui

	This file is part of the IGS software package for simulating
	the N-Body problem in real-time. This was developed as part
	of a departmental honors thesis project at Hofstra University
	in the Spring 2012 semester.

	IGS is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation; either version 2 of the License, or
	(at your option) any later version.

	IGS is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/

#include "Body.h"
#include "Gravity.h"
#include "NBodySystem.h"
#include <omp.h>

namespace IGS
{

#define SerialFor(i, min, max) for (int i = min; i < max; i++)
int MinParallelSize = 64;
int MinBlockSize = 32;

void BlockGravity(NBodySystem &bodies, const int i1, const int i2, const int j1, const int j2)
{
	const int ni = i2 - i1, nj = j2 - j1;
	if (ni > MinBlockSize && nj > MinBlockSize)
	{
		const int mi = i1 + nj / 2, mj = j1 + nj / 2;

		#pragma omp parallel sections
		{
			#pragma omp section
			BlockGravity(bodies, i1, mi, j1, mj);
			#pragma omp section
			BlockGravity(bodies, mi, i2, mj, j2);
		}

		#pragma omp parallel sections
		{
			#pragma omp section
			BlockGravity(bodies, i1, mi, mj, j2);
			#pragma omp section
			BlockGravity(bodies, mi, i2, j1, mj);
		}
	}

	else
	{
		SerialFor(i, i1, i2)
			SerialFor(j, j1, j2)
				bodies[i].Pairwise(bodies[j], bodies.Softening());
	}
}

void TriangleGravity(NBodySystem &bodies, const int i1, const int i2)
{
	SerialFor(i, i1, i2)
		SerialFor(j, i + 1, i2)
			bodies[i].Pairwise(bodies[j], bodies.Softening());
}

void BatchGravity(NBodySystem &bodies, const int n1, const int n2)
{
	const int n = n2 - n1;
	if (n > MinParallelSize)
	{
		int m = n1 + n / 2;

		#pragma omp parallel sections
		{
			#pragma omp section
			BatchGravity(bodies, n1, m);

			#pragma omp section
			BatchGravity(bodies, m, n2);
		}

		BlockGravity(bodies, n1, m, m, n2);
	}

	else
		TriangleGravity(bodies, n1, n2);
}

void SerialEvaluate(NBodySystem &system)
{
	TriangleGravity(system, 0, system.Size());
}

void ParallelEvaluate(NBodySystem &system)
{
	const int threads = omp_get_max_threads();

	// Allow batch parallelism when thread count is a power of 2
	if ((threads & (threads - 1)) == 0)
		BatchGravity(system, 0, system.Size());

	else
	{
		#pragma omp parallel for schedule(dynamic)
		SerialFor(i, 0, system.Size())
		{
			SerialFor(j, 0, system.Size())
			{
				if (i != j)
					system[j].Direct(system[i], system.Softening());
			}
		}
	}
}

}